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O Abstract 

(N 

Q>^ , We present a new algorithm for reducing an arbitrary unitary matrix U into a se- 

^ I quence of elementary operations (operations such as controlled-nots and qubit rota- 

tions). Such a sequence of operations can be used to manipulate an array of quantum 
bits (i.e., a quantum computer). Our algorithm applies recursively a mathemati- 
'^ ■ cal technique called the CS Decomposition to build a binary tree of matrices whose 

product, in some order, equals the original matrix U. We show that the Fast Fourier 



qh| Transform (FFT) algorithm is a special case of our algorithm. We report on a C++ 

program called "Qubiter" that implements the ideas of this paper. Qubiter(PATENT 
lj ■ PENDING) source code is publicly available. 



1. Introduction 

1(a) Previous Work 

In classical computation and digital electronics, one deals with sequences of 
elementary operations (operations such as AND, OR and NOT). These sequences 
are used to manipulate an array of classical bits. The operations are elementary in 
the sense that they act on only a few bits (1 or 2) at a time. Henceforth, we will 
sometimes refer to sequences as products and to operations as operators, matrices, 
instructions, steps or gates. Furthermore, we will abbreviate the phrase "sequence 
of elementary operations" by "SEO". In quantum computation[0], one also deals 
with SEOs (with operations such as controlled-nots and qubit rotations), but for 
manipulating quantum bits (qubits) instead of classical bits. Quantum SEOs are 
often represented graphically by qubit circuits. 

In quantum computation, one often knows the unitary operator U that de- 
scribes the evolution of an array of qubits. One must then find a way to reduce U 
into a SEO. In this paper, we present a new algorithm for accomplishing this task. 
We also report on a C-I--I- program called "Qubiter" that implements our algorithm. 
We call Qubiter a "quantum compiler" because, like a classical compiler, it produces 
a SEO for manipulating bits. Qubiter(PATENT PENDING) source code is publicly 
available at www.ar-tiste.com/qubiter.html. 

Our algorithm can be applied to any unitary operator U. 

It is useful to define certain unitary operators U^g for all A''^ G {1,2,3,...}, 
where U^g is a 2^^ x2^'^ matrix and Nb is the number of bits. Some U^g are known to 
be expressible as a SEO whose length (i.e., whose number of elementary operations) 
is a polynomial in Nb- Two examples are the Nb bit Hadamard Transform (HT) 
matrix and the Nb bit Discrete Fourier Transform (DFT) matrix. The HT matrix is 
known to be expressible as a SEO of length Order(A^B). The DFT matrix is known 
to be expressible (using the FFT algorithm of 0-0) as a SEO of length Order (A'^^). 
Our algorithm achieves both of these SEO-length benchmarks. Even better, the SEO 
often called the "quantum FFT algorithm" is exactly reproduced by our algorithm. 

We believe our algorithm yields short SEOs for many kinds of unitary operators 
other than the HT and DFT ones. We do not believe that our present algorithm 
yields the shortest possible SEO for every unitary operator. However, we do believe 
that it is possible to come very close to achieving this goal by introducing further 
optimizations into the algorithm. Future papers will report on our progress in finding 
such optimizations. 

Previous workers |^ have described another algorithm for reducing a unitary 
operator into a SEO. Like ours, their algorithm can be applied to any unitary operator 
U. However, it is very unlikely that their algorithm will be efficient in producing short 
SEOs unless further optimizations are added to it. And such optimizations, if they 
exist, have not been specified by anyone. Furthermore, as far as we know, there 



is no publicly available software that implements their algorithm. Our algorithm 
is significantly different from theirs. Theirs is based on a mathematical technique 
described in Refs.0-0, whereas ours is based on a mathematical technique called 
the CS Decomposition (CSD)[§-|jl2[ to be described later. 



Quantum Bayesian (QB) Nets|T^-[0] are a method of modeling quantum 
systems graphically in terms of network diagrams. In a companion paper [15 
show how to apply the results of this paper to QB nets. 



we 



1(b) CS Decomposition 

As mentioned earlier, our algorithm utilizes a mathematical technique called 
the CS Decomposition (CSD)|[-[[T2[. The C and S stand for "cosine" and "sine", 
respectively. Next we will state the special case of the CSD Theorem that arises in 
our algorithm. 

Suppose that U is a.n N x N unitary matrix, where A^ is an even number. 
Then the CSD Theorem states that one can always express U in the form 



U 



Lo 
Li 



D 



Ro 
Ri 



(Ib.l) 



where the left and right side matrices Lq,Li,Rq, Ri are "f- x y unitary matrices and 



D 






Doo = Dn = diag{Ci, C2, . . . , C, 



(lb.2a) 
(lb.2b) 



Doi = diag{Si, 6*2, . . . , Sjv) 



(lb.2c) 



Dio = -^01 ■ (lb.2d) 

For alH G {1, 2, . . . , ^1, Ci = cos 6i and Si = sin 9i for some angle 9i. Given any CSD 
of U, it is easy to find (see Appendix A) another CSD of U for which the angles 9i are 
in non-decreasing order and they are contained in the interval [0,90°]. Henceforth, 
we will assume that the angles 6i are so ordered and in this range. We will use the 
term D matrix to refer to any matrix that satisfies Eqs.(lb.2). If one partitions U 
into four blocks Uij of size ^ x y , then 



U^, 



^i ^ij ^j 7 



(lb.3) 

for i,j G {0, 1}. Thus, Dij gives the singular values [p!^] of Uij. 

More general versions of the CSD Theorem allow for the possibility that we 
partition U into 4 blocks of unequal size. 



Note that if U were a general (not necessarily unitary) matrix, then the four 
blocks Uij would be unrelated. Then to find the singular values of the four blocks Uij 
would require eight unitary matrices (two for each block), instead of the four Lj, Rj. 
This double use of the Li, Rj is a key property of the CSD. 

1(c) Bird's Eye View of Algorithm 

Our algorithm is described in detail in subsequent sections. Here we will only 
give a bird's eye view of it. 
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Fig.l A CSD binary tree. 



Consider Fig.l. We start with an initial unitary matrix Uin at horizontal level 
0. Without loss of generality, we can assume that the dimension of Uin is 2^^ for some 
Nb > 1. (If initially t/m's dimension is not a power of 2, we replace it by a direct 
sum Uin © Ir whose dimension is a power of two.) We apply the CSD method to Uin- 
This yields for level 1 a D matrix D{0, U), two unitary matrices L(0, U) and L(l, U) 
on the left side and two unitary matrices -R(0, U) and -R(l, U) on the right side. Then 
we apply the CSD method to each of the 4 matrices L{0,U),L{1,U),R{0,U) and 
R{1,U) that were produced in the previous step. Then we apply the CSD method 



to each of the 16 R and L matrices that were produced in the previous step. And 
so on. Each node M has children Afi and Mr to its left and right, respectively. Mr 
(ditto, Ml) has D, i?, L matrices produced by applying the CSD method to the R 
matrices (ditto, L matrices) of its parent node M . At level Nb-, the L's and i?'s are 
1x1 dimensional — i.e., just unit-modulus complex numbers. The nodes of the last 
level Nb + 1 don't have R,L or D matrices per se. If Mr (ditto. Ml) is a node of 
level Nb + 1, then it stores the R (ditto, L) matrices of its parent node M. 

Call a central matrix either (1) a single D matrix, or (2) a direct sum Di © 
-D2 © ■ ■ ■ © -Dr of D matrices, or (3) a diagonal unitary matrix. From Fig.l it is clear 
that the initial matrix Uin can be expressed as a product of central matrices, with 
each node of the tree providing one of the central matrices in the product. Later on 
we will present techniques for decomposing any central matrix into a SEO. 

2. Preliminaries 

In this section, we introduce some notation and some general mathematical 
concepts that will be used in subsequent sections. 

2(a) General Notation 

We define Za^i, = {a,a + 1, . . . ,b} for any integers a and b. 6{x, y) equals one 
ii X = y and zero otherwise. 

We will use the symbol Nb for the number (> 1) of bits and Ns = 2^^ 
for the number of states with Nb bits. Let Bool = {0, 1}. We will use lower case 
Latin letters a,b,c . . . G Bool to represent bit values and lower case Greek letters 
a, /?, 7, . . . e Zq^at^-i to represent bit positions. A vector such as a = a^^^i . . . 020109 
will represent a string of bit values, a^ being the value of the /z'th bit for /i G Zq^at^-i. 
A bit string a has a decimal representation d{a) = J2ij.=o^ '^^cl/i- For /? G Zqn^^i, we 
will use u{P) to denote the /?'th standard unit vector — i.e., the vector with bit value 
of 1 at bit position (3 and bit value of zero at all other bit positions. 

The set Bool^^ can be ordered in the standard way, the dictionary ordering. 
For example, the dictionary ordering of Bool^ is 

(000, 001, 010, Oil, 100, 101, 110, 111) . (2a.l) 

Other useful orderings of Bool^'^ are the so called Gray codes^^, named, not after 



the color, but after an actual person named Gray. In Gray codes, the next bit string 
may only differ from its predecessor in the value of a single bit. Clearly, this condition 
does not specify a unique ordering of Bool^^ , so there is more than one Gray code 
of Bool^'^ . Henceforth, we will refer to a Gray code as a lazy ordering., because, as 
we step from any b to the next, we act "lazily", flipping only one bit instead of many. 
An example of a lazy ordering of Bool^ is 



(000, 100, no, 010, on, in, loi, ooi) . 

We define the single-qubit states |0) and |1) by 



(2a.2) 



|0) 
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|1) 





1 



(2a.3) 



If a G Bool^'^ , we define the A^^-qubit state \a) as the following tensor product 



\a) = |aAr^_i) (g) . . . \ai) (g) \ao) 



(2a.4) 



For example. 
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(2a.5) 



Ir will represent the r dimensional unit matrix. Suppose /3 G Zq^tys-i and M 
is any 2x2 matrix. We define M{j3) by 



M(/5) 



M ®Io 



(2a.6) 



where the matrix M on the right side is located at bit position (3 in the tensor product 
of Nb 2x2 matrices. The numbers that label bit positions in the tensor product 
increase from right to left (^), and the rightmost bit is taken to be at position 0. 

For any pair of same-sized square matrices A and B such that det(-B) 7^ 0, one 
can define B^ = exp{A\n B) using the Taylor series of In(-) about the identity matrix 
and the Taylor series of exp(-) about zero. This gives a definition for objects such as 
Mi(/3i)^2(/32)^ where Mi and M2 are 2 x 2 matrices, det(Mi) ^ and (3i, p2 e ^o,7Vs-i- 

For any two same-sized square matrices A and -B, we define the o-dot product 
Qhy AQ B = ABA\ where A^ is the Hermitian conjugate of A. 

(^ = (c"a;, CTy, cT z) will represent the vector of Pauli matrices, where 



0-^ 



1 

1 



(T„ 



-i 
i 



cr. 



1 
-1 



(2a.7) 



A qubit rotation is defined as any matrix of the form exp[i^ ■ c?(/9)], where P G Zq^n^^i 
and ^ is a real 3-dimensional vector. 



2(b) Projection Operators 

Consider a single qubit first. 

The number operator n of the qubit is defined by 



n 




1 



1 - cr. 



Note that 



n|0) =0|0) = , n\l) = 1|1) 
We will often use n as shorthand for 



n = 1 — n 



1 




l + a. 



Define Pq and Pi by 



Pa = n 



1 




Pi = n 




1 



Pq and Pi are orthogonal projectors and they add to one: 

PaPb = S{ci, b)Pb for a,b E Bool 

P0 + Pl = l2. 



(2b.l) 



(2b.2) 



(2b.3) 



(2b.4) 



(2b.5) 



(2b.6) 



Now consider A''^ bits instead of just one. 

For /3 G Zq^Nj^_i, we define Po(/3), Pi{P), n{P) andn(/3) according to Eq.(2a.6).|T 

For a e Bool^'', let 



P^ = P 



a -^ ajsf I 



Pa2 ® Par ® Pao ■ 



For example, with 2 bits we have 

Poo = Po(^Po = dtag{l, 0,0,0) 

Poi = Po Pi = diag{0, 1, 0,0) 

Pio = Pi ® Po = diag{0, 0,1,0) 



P^^ = P-^^P-^ = diag{0, 0, 0, 1) 



Note that 



Nb 



P^P^ = S{d, b)P^ for a,be Bool 



^ p. = J2 ® J2 ® ■ ■ ■ (g) J2 = I2NB 
deBooi'^B 



(2b.7) 

(2b.8a) 
(2b.8b) 
(2b.8c) 
(2b.8d) 

(2b.9) 
(2b.l0) 
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For r > 1, suppose Pi, P2, • • • -Pr are orthogonal projection operators (i.e., 
PiPj = 6{i,j)Pj ), and ai, 02 • • • ctr are complex numbers. Then it is easy to show by 
Taylor expansion that 



exp(^ ttiPi) = J2 exp(ai)Pi + (1 - JI ^«' 



(2b.ll) 



In other words, one can "pull out" the summation sign from the argument of the 
exponential, but only if one adds a compensating term 1 — J2iPi so that both sides 
of the equation agree when all the ctj's are zero. 

For any pair of 2 x 2 matrices tq, ti, define 



T-b = rf, ®...0n,0no 



where b G Bool'^'^. Now let 



B{n,To) = {rAbeBool''^} 



(2b.l2) 



(2b.l3) 



The set of Ns x Ns diagonal complex matrices is a vector space and i3(Pi,Po) 
B{n,n) is a basis for it. B{az, h) and B{n, I2) are also bases for it. 



2(c) Sylvester-Hadamard Matrices 

The A^'^-bit Sylvester-Hadamard matrix |T^ Hnb is defined by: 
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(2c.la) 



(2c.lb) 



Hr+i = Hi®Hr, (2c. Ic) 

for any integer r > 1. We will often use a plain H to represent Hi. In Eqs.(2c.l), 
we have labelled the rows and columns with binary numbers in increasing dictionary 
order. From Eqs.(2c.l), one can show that the entry of Hn^ at row a G Bool'^^ and 
column b G Bool^^ is given by 



{H_ 



NB)d,b 



i-l 



^a-b 



(2c.2) 



where a ■ b = J2^=o '^^^/x- It is easy to check that 



Hnb = Hnb 1 (2c.3) 

Hl^ = NsIns ■ (2c.4) 

In other words, Hn^ is a symmetric matrix, and the inverse of H^g equals Hn^ 
divided by however many rows it has. 

If H{P) for P G Zq^at^-i is defined according to Eq.(2a.6), then Eqs.(2c.l) 
imply that 

H^g = H{Nb - 1) . . . H{2)H{1)H{0) . (2c.5) 

The if (/3)'s on the right side of the last equation commute so also 

H^g = H{0)H{1)H{2) . . . H{Nb - 1) . (2c.6) 

2(d) Discrete Fourier Transform 

The Nb bit Discrete Fourier Transform (DFT) matrix F/v^ is defined by 

iFNg)a,t = ^^"' , (2d.l) 

where a,b & Zq^Ns-i ci-nd 

u = e^ . (2d.2) 

Note that F^g is a symmetric matrix. It is easy to show that it is also a unitary 
matrix. If v and v are complex Ns dimensional vectors such that v = F^gV, then we 
call V the DFT of v. 

Calculating v the naive way, by multiplying v by -F/v^, would take 0rder(A^5^) 
classical elementary operations (complex multiplications mostly). Instead, it is pos- 
sible to calculate v from v in Order(A^s'ln A^s') classical elementary operations using 
the well known|0| Fast Fourier Transform (FFT) algorithm. Ref.[H was the first to 
express this algorithm as a product of matrices each of which acts on at most 2 bits 
at a time. This way of expressing it is often called the "quantum FFT algorithm" 
because it is ideal for quantum computation. Ref.p| showed (See Appendix B) that 

F^^ = -^H{Nb - 1) . . . A(2)ii(2) A(l)i7(l) A(0)ii(0)PsR , (2d.3) 

where H{a) is the 1-bit Hadamard matrix operating on bit a G Zo^at^-i, Pbr is the 
bit reversal matrix for A'"^ bits, and 

A(/3) = A(/5 + 1, /3) A(/3 + 2, /3) . . . A{Nb - 1, /3) , (2d.4) 



where 

A{a,p) = exp[i0|„_^|+in(a)n(/5)] , (2d.5) 

07 = ^ • (2d.6) 

Thus, A{a,f3) is a diagonal matrix whose diagonal entries are either 1 or a 
phase factor. For example, for Nb = 3, 

/n\ /oM \ J 1^2,01,00) if a2 = Oo = 1 MJ n'\ 

n(0)n(2)|a„ai,ao) = |Q ^^^^^^.^^ , (2d.7) 

so 

. , fn\ fo\ .000 001 010 Oil 100 ^^} 110 ^^} , 

A(0,2) = e*^3n(o)n(2) _ ^.^^^ 1,1,1,1, 1 ,e*^«, 1 ,e*'^^) . (2d.8) 

For Nb = 3, reversing the bits of the numbers contained in Zqj exchanges 
1 = rf(OOl) with 4 = rf(lOO) and 3 = rf(Oll) with 6 = rf(llO), and it leaves all other 
numbers in Zqj the same. Thus, for Nb = 3, Pbr is the 8x8 permutation matrix 
which corresponds to the following product of transpositions: (1,4)(3,6). 

Note that F^g, H{a), A{a) and Pbr are all symmetric matrices. Hence, 
taking the transpose of both sides of Eq.(2d.3), one gets 



F^g = -^PbrH{0)A{0)H{1)A{1)H{2)A{2) . . . H{Nb - 1) . (2d.9) 

Both the last equation and Eq.(2d.3) are called the quantum FFT algorithm. 
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• • 



H(N - 1) 



H(N„ - 2) 



H(N„ - 3) 



• • • 



H(N - 3) 



H(N - 2) 



H(N - 1) 



Fig. 2 Pictorial representation of quantum FFT algorithm. For simplicity, we only 
show the Hadamard matrices. White (unshaded) regions inside a matrix represent 

zero entries. 

Fig. 2 is a pictorial representation of Eqs.(2d.3) and (2d. 9). For any pair of 
2x2 complex matrices M,M' and for /?, /5' G Zq^Nb-i^ '^^ will say that M(/9) is 
"more diagonal" than M\j3) if /? < 13' . M(0) is a direct sum of 2 x 2 matrices. M(l) 
is a direct sum of 4 x 4 matrices so it is less diagonal than M(0), etc. The matrices 
A(a) that occur in Eqs.(2d.3) and (2d. 9) are truly diagonal but the matrices H[a) 
aren't. In Eq.(2d.3) (ditto, Eq.(2d.9)) the H{a) matrices become less diagonal (ditto, 
more diagonal) as one goes from right to left. 

The fact that the non-diagonal matrices H{a) are real and the diagonal matri- 
ces A(a) are complex is what gives the FFT algorithm its speed advantage in classical 
computation. Of course, segregating real and complex operations has no advantage 
in quantum computation, where complex numbers are very natural. The charm of 
the FFT algorithm in quantum computation is that it can be expressed as a short 
SEO, namely Eqs.(2d.3) and (2d.9). 
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2(e) Permutations 

Subsequent sections will use the following very basic facts about permutations. 
For more details, see, for example, Ref. [^ . 

A permutation is a 1-1 onto map from a finite set X onto itself. The set of 
permutations on set X is a group if group multiplication is taken to be function 
composition. Sn, the symmetric group in n letters, is defined as the group of all 
permutations on any set X with n elements. If X = Zi „, then a permutation G 
which maps i E X to ai E X (where i ^ j implies Oj 7^ aj) can be represented by a 
matrix with entries 

(G),, = 5(a„j), (2e.l) 

for all i,j G X. Note that all entries in any given row or column equal zero except 
for one entry which equals one. Hence, the rows of G are orthonormal and G is an 
orthogonal matrix {G^G = GG^ = 1). An alternative notation for G is 

G=f ^ 2 2 ■■■ " 1 . (2e.2) 

\ ai a2 a^ ■■ ■ an J 

The product of two symbols of the type shown in Eq.(2e.2) is defined by function 
composition. For example, 

I). (2e.3) 

Note how we have applied the permutations on the left side of the equation from 
right to left (<— ). (Careful: Some authors apply them in the opposite direction (—>)). 
A cycle is a special type of permutation, li G E Sn maps ai —^ a2, a2 —* a^, . . . , 
ar-i ^ Oj., Oj. — s> Oi, where i ^ j implies Oj 7^ a^ and r < n, then we call G a cycle. 
G may be represented as in Eqs.(2e.l) and (2e.2). Another way to represent it is by 

G = {ai,a2,a3,...,ar) . (2e.4) 

(Careful: some people write (a^, . . . , 03, a2, Oi) instead.) We say that the cycle of 
Eq.(2e.4) has length r. Cycles of length 1 are just the identity map. A cycle of length 
2 is called a transposition. The product of two cycles need not be another cycle. For 
example, 

(2,1,5)(1,4,5,6)=(^^ II 2 II) (2e-5) 

cannot be expressed as a single cycle. Any permutation can be written as a product 
of cycles. For example. 





12 3 4 5 6^ 
4 13 2 6 5^ 

12 



(5,6)(1,4,2). (2e.6) 



The cycles on the right side of Eq.(2e.6) are disjoint; i.e., they have no elements 
in common. Disjoint cycles commute. Any cycle can be expressed as a product of 
transpositions (assuming a group with > 2 elements), by using identities such as: 



[O-l, 0-2, 



an) = (ai, a2)(a2, as) ■ ■ ■ (a„_i, a„) , 



(2e.7) 



(ai, 02, ... , an) = (ai, a^) ■ ■ ■ (ai, a3)(ai, 02) . 
Another useful identity is 



(2e.8) 



{a,b) = {a,p){p,b){a,p) . (2e.9) 

This last identity can be applied repeatedly. For example, applied twice, it gives 



(a, 6) = (a,pi)(pi,6)(a,pi) = {a,pi){pi,p2){p2,b){pi,p2){a,pi] 



(2e.l0) 



Since any permutation equals a product of cycles, and each of those cycles can be 
expressed as a product of transpositions, all permutations can be expressed as a prod- 
uct of transpositions (assuming a group with > 2 elements). The decomposition of 
a permutation into transpositions is not unique. However, the number of transposi- 
tions whose product equals a given permutation is always either even or odd. An even 
(ditto, odd) permutation is defined as one which equals an even (ditto, odd) number 
of transpositions. 

3. State Permutations that Act on Two Bits 

The goal of this paper is to reduce any unitary matrix into a product of qubit 
rotations and controUed-nots (c-nots). A qubit rotation acts on a single qubit at a 
time. This section will discuss gates such as c-nots that are state permutations that 
act on two bits at a time. 



3(a) Nb = 2 

Consider first the case when there are only 2 bits. Then there are four possible 
states-00, 01, 10, 11. With these 4 states, one can build 6 distinct transpositions: 



(00,01) 
(00,10) 



C^x 




_ 


^2. 


'Pi 


Po' 


. ^0 


Pi _ 



PQ®a., + Pi®h = a,(0)"(i) 



(3a. la) 
(3a.lb) 
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(00,11) 



(3a. Ic) 



(01,10) 



(J-, 



(3a.ld) 



(01,11) 



(10,11) 



Po Pi 
Pi Po 



h (^ Po + (T^ (^ Pi = (Tx(l) 



cr^ 



n{0) 



Po®/2 + Pl®ax = ^^(0) 



n(l) 



(3a. le) 



(3a.lf) 



where matrix entries left blank should be interpreted as zero. The rows and columns 
of the above matrices are labelled by binary numbers in increasing dictionary order (as 
in Eq.(2c.lb) for H2). Note that the 4 transpositions Eqs.(3a.l)(a,b,e,f) change only 
one bit value. We will call them controlled nots (c-nots). The other 2 transpositions 
Eqs.(3a.l)(c,d) change both bit values. We will call (00, 11) the Twin-to-twin- er and 
(01, 10) the Exchanger. Expressions such as a^iP)"^""^ where a ^ (3 are a special case 
of Mi(/3i)*^2*^^2\ which was defined in Section 2a. ax{l3Y'^°'^ equals (Tx{I3) when it acts 
on a state for which n{a.) = 1, whereas it equals 1 if n{a) = 0. a is called the control 
bit and /5 the flipper bit. 

Exchanger [pl|] has four possible representations as a product of c-nots: 



(01, 10) = (01, 00)(00, 10)(01, 00) = a,(0)'^Wa,(l)"(°V,(0)"« 



(3a. 2a) 



(01, 10) = (10, 11)(11, 01)(10, 11) = a.(0)"«a.(l)"(°)a.(0)"« 



(3a.2b) 



(01,10) = (10,00)(00,01)(10,00) = a,(l)"(°V,(0)"«a,(l)"W 



(3a.2c) 



(01, 10) = (01, 11)(11, 10)(01, 11) = a.(l)"(°)a.(0)"«a.(l)"(°) . 



(3a.2d) 



Note that one can go from Eq.(3a.2a) to (3a. 2b) by exchanging n and n; from 
Eq.(3a.2a) to (3a. 2c) by exchanging bit positions and 1; from Eq.(3a.2a) to (3a. 2d) 
by doing both, exchanging n and n and exchanging bit positions and 1. We will 
often represent Exchanger by E{0, 1). It is easy to show that 
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E^(0,1) =^(0,1) = ^-^(0,1) , (3a.3a) 

E{0,1) = E{1,0) , (3a.3b) 

^2(0,1) = !. (3a.3c) 

Furthermore, if X and Y are two arbitrary 2x2 matrices, then, by using the matrix 
representation Eq.(3a.ld) of Exchanger, one can show that 

E(l,O)0(X®r) = r®X. (3a.4) 

Thus, Exchanger exchanges the position of matrices X and Y in the tensor product. 
Twin-to-twin-er also has 4 possible representations as a product of c-nots. One 
is 

(00,ll) = (00,01)(01,ll)(00,01) = (7,(Of«a,(ir(o)(T,(Of«. (3a.5) 

As with Exchanger, the other 3 representations are obtained by exchanging: (1) n 
and n, (2) bit positions and 1, (3) both. 
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A. 



ax(l)n(0)= >^ 



a^(l)n(0) = >^ 



n(l) = 



0^(0) 



1 1 

Fig. 3 Circuit symbols for the 4 different types of c-nots. 
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Fig. 4 Four equivalent circuit diagrams for Exchanger. 
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(00,11) = >^ 



y. 



i 
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y. 



= >^ 



y 
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A 



= 



t 



y 



X 



X = 



? 



o- 



Fig. 5 Four equivalent circuit diagrams for Twin-to-twin-er. 

Figures 3, 4 and 5 give a diagrammatic representation of the 6 possible trans- 
positions of states for Nb = 2. 

3(b) Any Nb>2 

Suppose ai, bi, 02, &2 £ Bool and a, P E Zq^at^-i such that a ^ 13. We define 



{aibi,a2b2] 



a,l3 



n 



(AaiA'foiA", AasA'fozA") 



(3b.l) 



(A,A',A")6-Boor 



where on the right side, 01,02 are located at bit position a, and 61, &2 are located 
at bit position f3. (Note that the transpositions on the right side of Eq.(3b.l) are 
disjoint so they commute.) For example, for Nb = 3, 

(j^(0)"(^) = (10,11)1,0= n («10'«11) = (010'011)(110'm) • (3b-2) 

aSiBool 

Clearly, any permutation of states with Nb bits that acts on only 2 bits (i.e.. Ex- 
changer, Twin-to-twin-er, and all c-nots) can be represented by (ai6i, a262)a,/3- 
For a,f3 E Zo^Af^-i such that a 7^ /?, let E{a,P) represent Exchanger: 

E{a,(3) = (01,10)^,0. (3b.3) 

As in the A^^ = 2 case, E{a, (3) can be expressed as a product of c-nots in 4 different 
ways. One way is 





y 



t 



^Wrr (RY<")rr ( r^V'W 



E(«,/?) = a,.(a)"^^V,(/3)"^"V,(a) 



(3b.4) 



The other 3 ways are obtained by exchanging: (1) n and n, (2) bit positions a and 
j3, (3) both. Again as in the Nb = 2 case. 



E^{a,f3) = E{a,(3) = E-\a,f3) 



(3b. 5a) 



E{a,P) = E{(3,a) 



(3b.5b) 
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E^{a,/3) = 1. (3b.5c) 

Furthermore, if X and Y are two arbitrary 2x2 matrices and a,/3 G Zq^^^^i such 
that a ^ (3, then 

E{a, (3) [X{a)Y{(3)] = X{(3)Y{a) . (3b.6) 

Equation (3b. 6) is an extremely useful result. It says that E{a,P) is a trans- 
position of bit positions. Furthermore, the E{a,P) generate the group of Nb\ per- 
mutations of bit positions. (Careful: this is not the same as the group of {2^'^)\ 
permutations of states with A''^ bits.) 



E(a, P) 



= i 



p a 



t T 

Fig. 6 Circuit symbol for Exchanger. 



X 



: 



>^ 



X 



X 



X 



▼ T t T T r 

Fig. 7 Circuit diagram for Eq.(3b.7). 
An example of how one can use Eq.(3b.6) is 



t 



E(2,0) [a,(0)"(i)a,(2)"(°V,(3)"(2)] = a,(2)"(iV,(0)"(2)a,(3)"(°) . 



(3b.7) 



Figure 6 gives a convenient way of representing E{a, (3) diagrammatically. Using this 
symbol, the example of Eq.(3b.7) can be represented by Fig. 7. 



18 



Certain experimental implementations of a quantum computer might only al- 
low nearest-neighbor bit interactions. Exchanger can be used to express any ele- 
mentary operation between 2 non-nearest-neighbor bits as a SEO that contains only 
nearest-neighbor interactions. For example, if A{3, 0) is an operator that acts on bits 
3 and 0, then 

A{3, 0) = [E{3, 2)E{2, 1)] A{1, 0) . (3b.8) 

Of course, identities that are true for a general transposition are also true for 
E{a,P). For example, 

(2,0) = (2,1)(1,0)(2,1). (3b.9) 

Therefore, 

E{2, 0) = E{2, l)E{l, 0)E{2, 1) . (3b.l0) 

4. Decomposing Central Matrix into SEO 

In Section 1(c), we gave only a partial description of our algorithm. In this 
section, we complete that description by showing how to decompose each of the 3 
possible kinds of central matrices into a SEO. 

4(a) When Central Matrix is a Single D Matrix 

D matrices are defined by Eqs.(lb.2). They can be expressed in terms of 
projection operators as follows: 

D= Y^ e-xp{i(t)d(yy) ® Pd , (4a.l) 

where the 0^ are real numbers. Note that in Eq.(4a.l), a has A^^ — 1 components 
instead of the full Nb- Using the identity Eq.(2b.ll), one gets 

(4a.2) 



D = exp 


i y ; (/)say (g) p. 




\ deBooi^B-'^ 


Now define new angles 6^ by 




0a = 


^ > ; i-ir% ■ 




baBooi^B-^ 



(4a.3) 

— * — * 

Suppose (ditto, 9) is a column vector whose components are the numbers 0^ (ditto, 
9^) arranged in order of increasing a. Then Eq.(4a.3) is equivalent to 

$=Hn^_,9. (4a.4) 
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This is easily inverted to 

— * 



0=7TW7ZtHns-i$- (4a.5) 



Let A^ for h G Bool^^ ^ be defined by 

Ag=exp(2%a,® Y. {-^f'Pa] ■ (4a.6) 

Then D can be written as 

D= n ^6- (4a. 7) 

b£Bool^B-^ 

Note that the A^ operators on the right side commute so the order in which they 
are muhiphed is irrelevant. Next we establish 2 useful identities: If /3 G Zo,Ar^_i and 
u{l3) G Bool^'^ is the /9'th standard unit vector, then 



j:,eBool-Bi-lf''^^^Pa 



^a&BoorB 

= /2 ® ■ ■ ■ ® /2 g 

= Po(/5)-Pi(/3) 
If /9, a G Zo,Afg_i and a ^ [3, then 



I]a^ eBooZ ( ~ 1 ) "'' -Pa^ 



(4a.8) 



a,(/3)"(")0a,(/3) 



a,(/5)[-Pi(a) + Po(a)] 
o-?;(/9)o-2(a) 



(4a.9) 



Now we are ready to express Ag in terms of elementary operators. For any h G 
Bool^^~^, we can write 

b = J2^iPj)^ (4a.l0) 

j=o 

where 

Nb -2> pr_i> ■■■> Pi > Po>0 ■ (4a.ll) 

In other words, b has bit value of 1 at bit positions f3j. At all other bit positions, b 
has bit value of 0. r is the number of bits in b whose value is 1. When b = 0, r equals 
0. By Eq.(4a.8) 

Ag = expli^ga,(iVB-l)n^.(/5j)) • (4a.l2) 
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By Eq.(4a.9), if r > 1, 

[a^{NB - l)"(^'-i) ■ --a^iNB - l)"(^^)a,(iVB - !)"(*)] ay{NB - 1) ., ^ „. 

Thus, for r > 0, 

Ag = [a^{NB - ir^f"-'^ ■ ■ ■ a^{NB - ir^f'^^a^iNj, - IT^^'^] exp[^%a,(iVB - 1)] , 

(4a.l4) 
where if r = 0, the expression to the left of is defined to be 1. There are other 
ways of decomposing Ag into a SEO. For example, using the above method, one can 
also show that 

Ag = [a^.^r^"^-'^ ■ ■ ■cT.(/3i)"(^^V,(/3o)"(^iV,(iV5 - 1)"(^«)] exp[2%a,(iVB - 1)] . 

(4a.l5) 
In conclusion, we have shown how to decompose a D matrix into a SEO. For 
example, suppose A''^ = 3. Then 



Define 6 by 



and r(-) by 



Then 



where 



D= Y. Q^v{i(t>ah(yy)®Pa®Pb- (4a.l6) 

a^baBool 



e = ^H2$, (4a. 17) 



r(^) = exp{i9(ry) 0/20/2. (4a. 18) 

D = AooAoiAioAn , (4a.l9) 

^00 = r(^oo) , (4a.20a) 

Aoi = a,(2)"(°)0r(^oi), (4a.20b) 

Aio = a,(2)"«0r(^io), (4a.20c) 

Au = [a,(2)"(iV,(2)"(°)] r(^n) • (4a.20d) 
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4(b)When Central Matrix is a Direct Sum of D Matrices 

Consider first the case A''^ = 3. Let R{(p) = exp(io"y0). Previously we used 
the fact that any D matrix D can be expressed as 

D= T. R{<P"a,)®Pa®Pb. (4b.l) 

a^baBool 

But what if R were located at bit positions or 1 instead of 2? The next two equations 
can be proven by expressing both sides of the equation as an 8 x 8 matrix. 

Do©A= E Pa®R{(t>'ab)®Pb, (4b.2) 

Dqo © Doi © Dio ®Dn= E Pa^Pb® R{(pab) , (4b.3) 

a,h£Bool 

where the Dj and Dij are D matrices. One can apply a string of Exchangers to move 
R in Eqs.(4b.2) and (4b. 3) to any bit position. Thus, 

Do®Di = E{1,2)q( J2 Wa6)®Pa®n) , (4b.4) 

\a,beBool I 



Doo © /^oi © ^10 ©Ai = [^(0,1)^(1,2)]© Y. R{(t>ab)®Pa®Pb\ . (4b.5) 

\a,b&Bool I 

(Careful: E(0,2) ^ E(0, 1)E(1, 2). E(0,2) will change Ea,fe^(0ab) ® Pa ® Pb to 
Yl,af> Pb® Pa® R{(pab), which is uot the same as the right side of Eq.(4b.3) ). 
For general Nb > 1, if /3 G Zq^n^^i and 



J 1 if /3 = or iVs = 1 

\ E{Nb -P-1,Nb-P)--- E{Nb -S,Nb- 2)E{Nb - 2, Nb - 1) otherwise 

(4b.6) 
then a direct sum of 2^ D matrices can be expressed as 

^0 ( E Ri(l>a)®Pa] . (4b.7) 

It follows that if we want to decompose a direct sum of D matrices into a SEO, 
we can do so in 2 steps: (1) decompose into a SEO the D matrix that one obtains 
by moving the qubit rotation to bit position Nb — 1, (2) Replace each bit name in 
the decomposition by its "alias" . By alias we mean the new name assigned by the bit 
permutation E defined by Eq.(4b.6). 



22 



4(c) When Central Matrix is a Diagonal Unitary Matrix 

Any diagonal unitary matrix A can be expressed as 

A= Y. exp(i03)P5, (4c.l) 

deBooi'^B 

where the 0^ are real numbers. Using the identity Eq.(2b.ll) yields 

A = expit J2 'PaPa] ■ (4c.2) 



Now define new angles 6^ by 



E (-l)'X- (4c.3) 

heBool^B 



In terms of vectors 



= HM^e , (4c.4) 

and 

^=^Hj,J. (4c.5) 

Let Ag for b G Bool^'^ be defined by 

\ aeSool^B I 

Then A can be written as 

A= n h^ (4c.7) 

h^Bool^B 

where the A^ operators commute. For any h G Bool^^ , we can write 

b = J2m)^ (4c.8) 

i=o 

where 

Nb-1> Pr-i > ■■■> Pi> Po>0 . (4c.9) 

(Careful: Compare this with Eq.(4a.ll). Now b G Bool^'^ instead of Bool^'^^^ and 
I3r-i can be as large as Nb — 1 instead of Nb — 2.) One can show using the techniques 
of Section 4(a) that 
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exp[i^o] if r = 
A^= { exp[i%cr^(/3o)] if r = 1 

[cT.iPor^'''-''^ ■ ■ ■ ^Mr^'^'^cT^iPor^"'^] © exp[z%a,(/3o)] if r > 2 

(4c.lO) 
As in Section 4(a), there are other ways of decomposing A^ into a SEO. 

In conclusion, we have shown how to decompose a diagonal unitary matrix 
into a SEO. For example, suppose Nb = 2. Then 

A = diagie''^'"', e'^'^^e'^^^ e^"^") . (4c. 11) 

Define 9 by 

d = -H2$. (4c.l2) 



By Eqs.(4c.7) and (4c. 10), 



where 



4 

A = AooAoiAioAn , (4C.13) 

Aoo = exp(i6'oo) , (4c. 14a) 

Am = h® exp(i^oiff^) , (4c. 14b) 

Aio = exp(i6'ioo-2) h , (4c. 14c) 

An = a.(0)"(') [h ® expiiOua,)] . (4c.l4d) 



4(d) Comments 



The use of a Hadamard transform in Section 4(a) may seem at first somewhat 
mysterious to the reader. Sorry. Here is some motivation for it. 

We began Section 4(a) by noting with Eq.(4a.2) that any D matrix D can be 
expressed as -B = exp[iay{NB — 1)] raised to some power which equals a linear com- 
bination of products of PoiP) = n{P) and Pi{(3) = n{P), where P G Zo,Afs-2- Using 
a Hadamard transform allowed us to re-express D as a product of Ag's . According 
to Eq.(4a.l2), each A^ can be expressed as B raised to some power which equals a 
product of (Jz{P), where P G .^o,Afs-2- If in Eq.(4a.2) we had replaced each n(/9) and 
n{P) by |[1 — o-ziP)] and |[1 + <7z{P)], respectively, then we would have found the 
same thing that we found via a Hadamard transform — i.e., that D is expressible as a 
product of ^g's, each of which is expressible as B raised to some product of aziPYs. 
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For any 2x2 complex matrix B such that det(-B) 7^ 0, we define a controlled 
gate g to be an operator of the form 

g = B{af'o(MPH(M-P.^.,(A^^) ^ (4^.1) 

where bo, &i, • • • , K-i G Bool, and Po, Pi, ■ ■ ■ , Pr-i, a are distinct elements of Zq^n^^i. 
We call Po, Pi, . . . , Pr-i the control hits and a the flipper hit. An example of Eq.(4d.l) 
is when B = e*"^ for some real 0, in which case we call g a controlled phase factor. We 
already encountered controlled phase factors, with 2 controls, in the FFT algorithm. 
Two other examples of Eq.(4d.l) are when B = crx{<y) and B = exp[z6' ■ 0^(0;)] for some 
real 3-dimensional vector 6, in which cases we get a c-not (with possibly more than 
one control) and a controlled quhit rotation. 

In general, whenever one has a controlled gate g with more than one control 
bit, one can simplify g by expressing its exponents Po{P) = n{P) and Pi{P) = n{P) in 
terms of hiP) and (Jz{P). Appendices C and D say more about this transformation. 
Appendix D discusses the transformation from the point of view of Linear Algebra. 
Appendix C uses the transformation to decompose a controlled gate with 2 controls 
into a SEO. Decompositions like the one in Appendix C were first discovered by the 
authors of Ref.[^. However, their method of deriving such decompositions is very 
different from ours. 



5. SEO-Length Optimizations 

In previous sections, we have given an algorithm that can decompose an ar- 
bitrary unitary matrix into a SEO. But the SEO's generated by this algorithm may 
not be the shortest possible. In this chapter we will indicate several "SEO-length 
optimizations" — small adjustments to the algorithm which are either guaranteed or, 
at least, likely to produce shorter SEO's. 

The group U{Ns) of A^^ x Ns unitary matrices has Ns"^ free (real) parameters. 
(It has Ns"^ complex entries for a total of 2Ns'^ real parameters, but those parameters 
must satisfy Ns'^ orthonormality constraints.) In a CSD tree, level r has 2^~^ nodes 
so the number of nodes in the tree is 1 + 2 + 2^ + . . . + 2^^ = 2^^+! _ 1 ~ 2Ns. (Level 
Nb + 1 alone has Ns of the nodes, about half of them!) Each node yields Order (A^g) 
elementary operations. Therefore, in general, a CSD tree yields a SEO whose length 
is Order(A^s'^). In essence, each elementary operation of the SEO (except permutation 
operations like c-nots) carries (approximately) one of the Ns'^ free parameters. 

But we know that certain families of matrices contained in U{Ns) have internal 
symmetries which allow us to parametrize them with substantially fewer than Ns"^ 
parameters. And we know how to decompose some of these families (for example, the 
DFT matrices) into a SEO whose length is a polynomial in Nb. Can such SEOs be 
obtained with our CSD tree algorithm? As we shall see, the answer is yes, at least in 
some important cases like the DFT matrices. The trick is to introduce optimizations 
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into our algorithm which make the CSD tree degenerate into a simple string (or 
nearly one) of nodes. The other nodes never "grow" . For any node A/", if its right side 
matrices(i.e., i?(6,A/') in Fig.l, where b G Bool^ and A is the level of A/") are all equal 
to the identity matrix, then M bears no children on its right side. Ditto for its left 
side matrices. Note that in Fig. 2, the first (ditto, second) quantum FFT expansion 
looks like what one would expect if a CSD tree were to degenerate into the string 
consisting of the rightmost (ditto, leftmost) node of each level of the tree. As we shall 
see, this is precisely what happens if we add the following optimizations to our CSD 
tree algorithm. 

5 (a) Lazy Ordering of Factors 

Consider Eq.(4a.l9). The operators A^ commute so this equation is valid 
regardless of the order in which the A^ are multiplied. Suppose that we multiply the 
A^ so that their subscripts h are in a lazy ordering that starts with 00 on the right 
and ends with 01 on the left: 

D = AoiAnAioAoo . (5a. 1) 

After inserting Eqs.(4a.20) for the Ag, this yields 

D = a.(2)"(°)r(^oi)^.(2)"(^)r(^n)^.(2)"(°)r(^io)^x(2)"(^)r(^oo) . (5a.2) 

We see that by ordering the Ag in this way, several c-nots cancel out. Only one c-not 
remains between adjacent F's. There is no c-not to the right of the rightmost F in 
Eq.(5a.2), because we started the lazy ordering with b = 00. There is only one c-not 
to the left of the leftmost F in Eq.(5a.2), because we ended the lazy ordering with 

6 = 01. 

The above example assumes Nb = 3 and that we are decomposing a central 
matrix of type 1 (i.e., a central matrix which is a single D matrix). However, this 
method of shortening a SEO can also be used for the other two types of central ma- 
trices and for other values of Nb- Indeed, the other two types of central matrices are 
also decomposed by the algorithm into a product of Ag type matrices. Furthermore, 
it is possible to find for arbitrary Nb a lazy ordering of Bool^^ which starts with the 
zero vector 6 = and ends with a h which has only one non-zero component. 

5(b)Lightening Right Side Matrices 

Suppose D is a 4 X 4 D matrix such that its two angles are equal. Hence, 



D 
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(5b.l) 



where c = cos^, s = sin6 for some angle 6. Then for any 2x2 unitary matrix Gq, 

(Go © Go)D{Gl Q)Gl) = D . (5b.2) 

More generally, consider an A^ x A^ D matrix D. Let 6a for a ^ Z^ n__^ be its 
angles, arranged in non-decreasing order (recall Appendix A). Suppose the first do 
angles are equal to each other, the next di angles are equal to each other but larger 
than the previous di angles, and so on. Hence, do + di + . . . + dpi-i = y- We will 
call do,di, . . . , du-i the degeneracies of the D matrix angles. For a G Zo^m-I) let Ga 
be any da x da unitary matrix, and define G by 



Then 



G = Go®Gi® ...® Gm~i ■ (5b.3) 



{G ® G)D{G^ ® G^) = D . (5b.4) 



Suppose that 



U = {Lo®Li)D{Ro®Ri) (5b.5) 

is a CSD of a unitary matrix U . If G is a unitary matrix that satisfies Eq.(5b.4), 
define L'^ and B!^ for j G Bool by 

L] = LjG , (5b.6) 

R'^ = G^Rj . (5b.7) 

U = (L'o © L[)D{R'q © R[) (5b.8) 



Then 



is also a CSD of U. Suppose D has M distinct angles with degeneracies do,di, . . . du-i- 
Define S_i = and T^a = do + di + ■■■ + da ioi all a G Zq^m-i- For a G Zq^m-i and 
j G Bool, suppose R!j (ditto, Rj ) are the da rows with row indices from T^a-i to 
Sa — 1 of R!j (ditto, Rj). Then Eq.(5b.7) is equivalent to the equations 

Rf = GaR'l"^ , (5b.9) 

for a G Zo^M-i and j G Bool. For each a G Zq^m-i, "we will choose the two matrices Ga 
and R'o on the right side of the last equation so that they are the Q and R matrices. 



respectively, in a QR decomposition of Ro ■ In a QR decomposition [jlGj , the Q matrix 



is unitary so constructing Ga this way will make it unitary, as required by our original 
definition of it. For any rectangular matrix M, its principal diagonal are those entries 
of M for which the row and column indices are equal. In a QR decomposition, the R 
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matrix has zero entries below the principal diagonal so Rq will have this property. 
In a QR decomposition, it is possible to multiply each row of the R matrix by a 
phase factor and the corresponding column of Q by the conjugate phase factor. Thus 
it is possible for each a G Zqm-i to do a QR decomposition in which all the terms 
along the principal diagonal of Rq are non-negative. (Alternatively, one can make 
non- negative all the terms along the principal diagonal of -Rq.) 

By lightening the right side matrices, we will mean replacing the CSD Eq.(5b.5) 
by the CSD Eq.(5b.8), where G is chosen so that for all a E Zq^m~i, the Rq have 
zero entries below the principal diagonal and they (or else Rq) all have non-negative 
diagonal entries. The strategy of lightening the right side matrices is not guaranteed 
to make a CSD tree have fewer nodes. However, whenever it does achieve this goal, 
the reason why it does so is clear. The strategy tries to make zero as many entries 
of Rq as possible. This causes the CSD tree to grow mostly to the left. Its right side 
growth is stunted. 

5(c) Extracting Phases From Complex D Matrices 

For this optimization, it is convenient to generalize the definitions of a D 
matrix and of a central matrix. Henceforth, let a D matrix be any unitary matrix D 
which can be written in the form 



D 



Dqo Dqi 



(5c. la) 



where the Dij for i,j G Bool are same-sized diagonal matrices. An equivalent defini- 
tion: a D matrix is a matrix of the form 

D= Y. f^^"^®^a, (5c. lb) 

deBooi^B-'^ 

where the U^"'' are unitary 2x2 matrices. What we previously called a D matrix 
will now be called a real D matrix. If we need to emphasize that a D matrix is of the 
new kind, we will describe it as a complex D matrix. Henceforth, let a central matrix 
be defined as before, except that the term "D matrix" in the previous definition of 
a central matrix should now be interpreted to mean a complex D matrix. What we 
previously called a central matrix will now be called a bare central matrix{we can't 
call it a real central matrix because it could be a complex unitary diagonal matrix). 
If we need to emphasize that a central matrix is of the new kind, we will describe it 
as a dressed central matrix. 

In previous sections, we showed that: (1) Given an arbitrary unitary matrix 
U, one can construct a CSD tree in which each node stores a bare central matrix 
such that U equals the product (in a particular order) of these central matrices. (2) 
Any bare central matrix can be decomposed into a SEO. In this section, we will show 
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(1) and (2), but with the words "bare central matrix" replaced by "dressed central 
matrix" . 

To construct a CSD tree of dressed central matrices, one proceeds as follows. 
Suppose node M at level A G Zi^tvs-i has children nodes Ml and Mr to its left and 
right, respectively. Node M will store a dressed central matrix Q{M), a collection 
of left side matrices L{b,M) for b G Bool^, and a collection of right side matrices 
R{b,M) for b G Bool^. One performs a normal CSD only on those side matrices 
R{b,M) and L{b,M) which are not complex D matrices. If a side matrix is a complex 
D matrix, one performs an aborted CSD. Given a subroutine f{U,LQ,Li,D,RQ,Ri) 
for doing CSD, with input U and outputs Lq, Li, D, Rq, Ri, an aborted CSD is defined 
as a pass through the subroutine that returns the identity matrix for the side matrices 
Lq, Li, Rq, Ri and it returns U for D. 

So far we have shown (1) — how to construct a CSD tree of dressed central 
matrices. Next we will show (2) — how to decompose a dressed central matrix into a 
SEO. We will show how any dressed central matrix Q can be factored into 

Q = AlQ'Ar , (5c.2) 

where Q' is a bare central matrix with real entries and where A^ and Ar are diagonal 
unitary matrices (so they are also bare central matrices, though not necessarily real). 
In previous sections, we showed how to decompose any bare central matrix into a 
SEO. So we will be done proving (2) if we can show Eq.(5c.2) 

Consider a 2 x 2 complex D matrix D. One can always express D as 



D 






(5c.3) 



where fl, ul, ojr, c, s are real numbers, and where c^ + s^ = 1. Eq.(5c.3) implies 



D = diag{\, e 



lUJl^\ 



c s 
—s c 



e'''diag{l,e'^'') . (5c.4) 



Note that we can always choose s > 0. For if s < 0, then one can replace s -^ \s\, 
ujl ^ ujl + 71 and ujr -^ ur + vr. Likewise, we can always choose c > 0. For if c < 0, 
then one can replace c — ^ |c|, cj^ ^ cj^ + vr, ur -^ uor + tx and Vt ^ Vt + -k. With 
s > and c > 0, we can write c = cos^ and s = sin 6 with 6 G [0, 90°]. It is a trivial 
exercise to solve for the parameters Q, ujl, ojr and 9 in terms of the entries of D. 
Now consider an N x N complex D matrix D. Let 

F,, = rfm^(e*^^"'e*^^'\ . . . , e''^^^''^) for a G {L, R} , (5c.5) 

F = rfm^(e^^'°^e^^^^...,e^^*^"V (5c.6) 
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C = diag{Co, Ci, . . . , Cn_^) , (5c. 7) 

S = diag{So, Si, ... , Sn_^) , (5c. 8) 

where for any j E Zq n__^ and any a E {L, R}, fi*^-'-' and cu^^ are real, and Cj = cos 9j, 
Sj = sinOj with 6j E [0,90*^]. One can always express D as 



D = [In® Tl] 



c s 
-s c 



[r © (rr^)] . (5c.9) 



Note that this is really a special CSD. To harmonize with our similar CSD angle 
convention (see Appendix A), it is important that we choose the angles 6j to be in 
the interval [0,90°]. As in the N = 2 case, it is trivial to solve for the parameters 
Q^U)^ (jU) and 9j in terms of the entries of D. 

Now consider any dressed central matrix Q. We want to show that Q = 
AlQ'Aji. This equation is trivially satisfied if Q is a diagonal unitary matrix so 
assume instead that Q is a direct sum of one or more complex D matrices. For 
example, suppose that 

Q = D{00) © D{01) © D{10) © D{11) , (5c.l0) 

where the matrices D{b) for b E BooP are complex D matrices. By Eq.(5c.9), for 
each b E Bool"^, one can write 

D{b) = AL{b)D'{b)AR{b), (5c.ll) 

where the matrices Ai(6) and Aji{b) are diagonal unitary matrices and D'{b) is a real 
D matrix. Q = A^Q' Ar now follows if we set 

A, = A,(00) © A,(01) © A,(10) © A,(ll) , (5c.l2) 

for a E {L, R} and 

Q' = D'{00) © D'{01) © D'{10) © D'{11) . (5c.l3) 

One can expand A^, Q' and Aji into a SEO using the techniques of Section 
4. Alternatively, one might do this only for Q', and express A^ and Ar as a product 
of controlled phase factors. We end this section with a discussion of controlled phase 
factors and of their products. 

Any Ns X Ns diagonal unitary matrix A can be expressed as a product of 
controlled phase factors. For example, suppose 

A = diag (e^'^«» , e''^"' , e^'^^" , e^"^" ) , (5c. 14) 
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where 0g is a real number for all b G Bool^. Eq.(5c.l4) can be rewritten as 

By replacing all exponents n(/3) with 1 — n{[3), the last equation can be expressed in 
the form: 

where the 6^ are real numbers, linear combinations of the 0g. The SEO of Eq.(5c.l6) 
is preferable to that of Eq.(5c.l5) because Eq.(5c.l5) contains four 2-bit operations 
whereas Eq.(5c.l6) contains only one. From the point of view of Linear Algebra, 
going from Eq.(5c.l5) to Eq.(5c.l6) is an example of changing from the B{n,n) to 
the B{n, 1) basis. See Appendix E for more information about this transformation. 

In this paper, we use bare central matrices whose D matrices have angles 
contained in the interval [0, 90°]. An interesting special case is when these angles are 
all either or 90° exclusively. The entries of such central matrices are elements of 
{—1, 0, 1}. When the optimization being discussed in this section is turned ON, it is 
sometimes convenient to decompose such central matrices into a SEO by a method 
different from the one discussed in Section 4. See Appendix F for more information 
about this. 

5(d) Permuting Bits before Each CSD 

This optimization is only partially implemented in the current version of 
Qubiter, and it was not used to get any of the results of Section 6. 

level Ui„ (not a node) 



level 1 



level 2 




level 3 ^ // = ULR 

Fig.8 Subset of CSD tree. 
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We will start with an example and then generalize it. Suppose that A^^ > 3 
and that we are interested in finding the central matrix belonging to node A/3 = ULR 
of Fig. 8. This being a binary tree, A/3 has only one parent, call it A/2 = UL. A/2 in 
turn has just one parent A/i = U, which is the level 1 node (root node) of Fig.8. If Um 
is the initial unitary matrix and P(0,A/i) is some permutation matrix that permutes 
bits (hence, P(0, A/i) is always expressible as a product of Exchangers), we do a CSD 
of P{0,J\fi) Uin instead of ?/,„: 



P(0, ATi) U^n = [L{0,Mi) © L{l,^i)]D{0,M,)[R{0,^i) © i?(l, A^i)] . (5d.la) 

To get the central matrix belonging to A/2, we permute the bits of L(0,A/i) and 
L(l, A/i), and then do a CSD of the outcome: 



P{a,M2) L{a,Mi) = [L{aO,M2) © L{al,M2)]D{a,M2)[R{aO,M2) © R{al,^2)] , 

(5d.lb) 
for all a G Bool. Finally, to get the central matrix belonging to A/3, we permute the 
bits of -R(a6, A/2) for all a,b & Bool, and then do a CSD of the outcome: 



P{ah,M^)QR{ah,M2) = [L{abO,^f^)®L{abl,^f3)]D{ab,^f3)[RiabO,^f3)®Riabl,^fs)] , 

(5d.lc) 
for all a,b E Bool. The central matrix Q{N'j) for j G Zi^ is given by 

Q{Ui) = D{Q,Mi) , (5d.2a) 

Q{N2) = DiO,X2) © Dil,X2) , (5d.2b) 

QiAfs) = 0(00,^3) ® 0(01,^3) ® 0(10,^3) ® D{11, Ms) . (5d.2c) 

One can define a permutation matrix P{Mj) for j G ^1,3 by 

P(Ari) = P(0,Ari) , (5d.3a) 

P{M2) = P(0, Ar2) © P(l, Ar2) , (5d.3b) 

P(Ar3) = P(00, ATg) © P(01, ATs) © P(10, Ar3) © P(ll, Ar3) . (5d.3c) 

It is useful to define matrices Q'{A/j) for j G ^1,3 by undoing the permutations that 
were performed to calculate Q'{Afj): 
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Q'(M) = P-\^fl) Q(M) , (5d.4a) 

Q\Af2) = [PiAf2)P{Afi)]-' Q{Af2) , (5d.4b) 

Q'i^f3) = [Pi^fs)Pi^f2)Pi^fl)]-' giM) . (5d.4c) 

When the optimization presently being discussed is turned OFF, we associate 
with each node A/" a central matrix Q{J\f), and the product of these Q{M) in some 
special order equals f/j„. With this optimization turned ON, we associate with each M 
a pair {Q{M)^ P{N')) consisting of a central matrix Q{M) and a permutation matrix 
P(A/'), and the product of the Q'i^hf) equals f/j„. Note that the matrices Q'{Af) are 
not necessarily central matrices. 

The above discussion can be generalized to arbitrary A''^ as follows. Assume A/a 
is a node at level A G -Z'i^tvs- A/a has a string of predecessors A/a_i, A/a_2, • • • A/i, where 
^fj-l is the single parent oi Afj and A/i is the tree's root node. Define BooP = {0}. 
If A = 1, let Ux-i = Uin- If A > 1, let J/a-i equal either R(b,J\fx-i) or L(6,A/'a-i) 
for some b G Bool^^^. For any 6 G Bool^^^, the D,R,L matrices of node A/a are 
obtained by applying to Ux^i a permutation matrix P{b,Af\) that permutes its bits 
and then doing a CSD on the outcome: 

P(6, ATa) f/A-i = [L(60, ATa) © L(61, ArA)]D(6, AAa)[P(W, ATa) © i?(61, AT,)] . (5d.5) 
Define the central matrix Q{Afx) of node A/a by 

Q(A/'a)= D{b,^fx). (5d.6) 

Define the permutation matrix P(A/a) of node A/a by 

nA/'A)= P(fe,ArA). (5d.7) 

Finally, define the matrix Q^Mx) of node A/a by 

Q\Ux) = [P{Ux) ■ ■ ■ P{U2)P{Ui)r^ Q{Ux) ■ (5d.8) 

With this optimization ON, we associate a pair {Qi^M), P{Af)) with each node 
A/" of the CSD tree, and the product of the (not necessarily central) matrices Q'{J^) 
in some order equals f/j„- 

How to best choose the permutation matrices P(A/') for this optimization re- 
mains an open question. One can carry out two main types of permutation searches: 
exhaustive or heuristic. An example of a partially exhaustive permutation search 
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(and the only permutation search implemented in the current version of Qubiter) is 
as follows. Try all possible bit permutations -P(A/i) for the root node. Set -P(A/') equal 
to the identity for all other nodes. Of all possibilities tried, find that which yields the 
shortest SEO. Obviously, one can try more complicated kinds of exhaustive searches. 
As for heuristic searches, for each node M at level A, one could try to find a permu- 
tation -P(A/') which produces the smallest possible number of distinct CSD angles in 
the matrices D{b,f/) for all b G Bool^~^ . One expects that increasing the degeneracy 
in these CSD angles will make the "lightening right-side matrices" optimization more 
effective. 

6. Qubiter 

At present, Qubiter is a very rudimentary program. The current version — 
Qubiterl.l — is written in pure C++, and has no graphical user interface. All opti- 
mizations discussed in Section 5 except 5(d) have been fully implemented in Qubiterl.l. 
In its "compiling" mode, Qubiter takes as input a file with the entries of a unitary 
matrix and returns as output a file with a SEO. In its "decompiling" mode, it does 
the opposite: it takes a SEO file and returns the entries of a matrix. The lines in a 
SEO file are of 6 types. 

(a) ROTY a ang 

where a G Zq^Nb-i ^''^^ '^^fi' is a real number. This signifies the rotation of 
qubit a about the Y axis by an angle of ang in degrees. In other words, 
exp{iay{a)ang^). (Some people would call this a rotation by 2 ang instead 
of ang). 

(b) ROTZ a ang 

This is the same as (a) except that the rotation is about the Z axis instead of 
the Y one. 

(c) SIGX a 

where a G Zq^Nb-i- This signifies unconditional fiipping (NOT) of qubit a. In 
other words, ^"^^(a;). 

(d) CNOT «! chari a2 char2 ... ftr chavr [3 

where r > 1, ai, 0:2, • • • «r, /5 are distinct elements of Zq^^b-i ^ind chari, char2, ■ ■ ■ chavr 
are elements of {T, F}. This signifies a controlled- not with r controls. First sup- 
pose r = 1. If chari is the character T, this signifies aj.{P)"'^°''^\ Read it as: 
"c-not: if ai is true, then fiip (3." If chari is the character F, this signifies 
o'xiP)"'^"^^ ■ Read it as: "c-not: if ai is false, then fiip /3." Cases with r > 1 are 
defined analogously. For example, CNOT T 1 F 2 signifies ax{2)^^^^''^^^\ Read 
it as: "c-not: if bit is true and bit 1 is false, then fiip bit 2." 
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(e) PHAS ang 

where ang is a real number. This signifies a phase factor exp(i{ang)^). Thus, 
ang is an angle expressed in degrees. 

(f ) CPHA ai chari a^ char2 ... Ur chavr ang 

where r > 1, ai, 02, . . .ar are distinct elements of Z^at^-i, chari, char2, . . . c/iar^ 
are elements of {T,F}, and ara^f is a real number. This signifies a controlled 
phase factor with r controls. For example, PHAS T 1 F 90 would signify 
[exp(i90^)]"(°)"(i). Read it as: "phase shift: if bit is true and bit 1 is false, 
then shift phase by 90 degrees." Thus, ang is an angle expressed in degrees. 

CNOT with more than one control and CPHA with more than two are not elementary; 
i.e., they act on more than 2 bits. Strictly speaking, non-elementary operations should 
not be allowed in a SEO. However, Qubiter allows these 2 types because they are 
popular in the literature. We have indicated in previous sections how to decompose 
such controlled gates into SEOs. 

Fig. 9 shows the output of Qubiter for the 2, 3 and 4 bit Hadamard matrices, 
with all optimizations mention in Section 5 except 5(d) turned ON. Note that Qubiter 
is smart enough to realize that it is dealing here with a tensor product of 1 bit 
matrices. Fig. 10 shows the output of Qubiter for the 2, 3 and 4 bit DFT matrices, 
with all optimizations mentioned in Section 5 except 5(d) turned ON. Note that the 
quantum FFT algorithm of Ref.[^ is exactly reproduced. Thus, the quantum FFT 
algorithm is a special case of our CSD tree algorithm. In both the Hadamard and 
DFT cases, the CSD tree degenerates into the string consisting of the leftmost node 
of each level of the tree. 

In the future, we plan to introduce into Qubiter more optimizations and some 
quantum error correction. Much work remains to be done. 
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2 BITS: 

ROTY 1 45.0000000 
CPHA 1 T 180.000000 
ROTY 45.0000000 
CPHA T 180.000000 



3 BITS: 

ROTY 2 45.0000000 

ROTY 1 45.0000000 

CPHA 1 T 180.000000 

CPHA 2 T 180.000000 

ROTY 45.0000000 

CPHA T 180.000000 



4 BITS 


ROTY 


3 


ROTY 


2 


ROTY 


1 


CPHA 


1 


CPHA 


2 


CPHA 


3 


ROTY 





CPHA 






45.0000000 
45.0000000 
45.0000000 
T 180.000000 
T 180.000000 
T 180.000000 
45.0000000 
T 180.000000 

Fig. 9 Output of Qubiter with input matrix equal to -^11^^ for Nb = 2, 3, 4. H^^ 
is the Nb bit Hadamard matrix defined in Section 2(c). 
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2 BITS 




ROTY 


1 


45.0000000 


CPHA 


1 


T 180.000000 


CPHA 


1 


TOT 90.0000000 


ROTY 





45.0000000 


CPHA 





T 180.000000 


3 BITS; 




ROTY 


2 


45.0000000 


CPHA 


2 


TIT 90.0000000 


ROTY 


1 


45.0000000 


CPHA 


1 


T 180.000000 


CPHA 


1 


TOT 90.0000000 


CPHA 


2 


T 180.000000 


CPHA 


2 


TOT 45.0000000 


ROTY 





45.0000000 


CPHA 





T 180.000000 


4 BITS; 




ROTY 


3 


45.0000000 


CPHA 


3 


T 2 T 90.0000000 


ROTY 


2 


45.0000000 


CPHA 


2 


TIT 90.0000000 


CPHA 


3 


TIT 45.0000000 


ROTY 


1 


45.0000000 


CPHA 


1 


T 180.000000 


CPHA 


1 


TOT 90.0000000 


CPHA 


2 


T 180.000000 


CPHA 


2 


TOT 45.0000000 


CPHA 


3 


T 180.000000 


CPHA 


3 


TOT 22.5000000 


ROTY 





45.0000000 


CPHA 





T 180.000000 



Fig. 10 Output of Qubiter with input matrix equal to PbrFnb for ^b = 2, 3, 4. F^g 
is the Nb bit DFT matrix defined in Section 2(d). 

Appendix A. We Can Assume CSD Angles Are 
Non-decreasing Elements of [0, 90^ 



In this appendix we will show that: given any CSD of U as in Eq.(lb.l), one 
can always find another CSD of U for which the angles 9i are in non- decreasing order 
and they are contained in the interval [0, 90°]. 
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For simplicity, suppose the number of CSD angles is two. The discussion 
that follows can be easily generalized to any number of angles. For two CSD angles, 
Eq.(lb.l) becomes 



U 



U 



where 
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Si 


" 






"i^o " 





C2 


S2 


"i?o 





Li 


-si 


Cl 








Ri 







-S2 
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0. Then 










Cl 


si 


" 






"^0 1 





|C2| 


-S2 


' R, 





l; 


-si 


Cl 








Ri 







S2 


C2 







(A.l) 



(A.2a) 



L'- = Ljdiag{l, —1] 

for j = 0, 1. Hence, we may assume ci, C2 > 0. 
Suppose Cl, C2 > but S2 < 0. Then 



(A.2b) 



where 



U 









Cl 





Sl 


" 








Lo 










C2 





S2I 




' Ro 








L[_ 




-si 





Cl 










R'l 











- S2I 





C2 . 









(A.3a) 



L[ = Lidiag{l, —1] 



(A.3b) 



R[ = diag{l, —l)Ri 



Hence we may assume ci, C2, si, S2 > 0. 

Finally, suppose that 61 > 62- Then 



U 



where 



L'o 
L[ 



C2 





S2 








Cl 





Sl 


-S2 





C2 








-Sl 





Cl 



i?'o 
R[ 



(A.3c) 



(A.4a) 



L', = Li 



1 

1 



(A.4b) 
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i?;= J J i?,, (A.4c) 

for j = 0, 1. Hence, we may assume that 6i < 62- 

Appendix B. Proof of Quantum FFT Algorithm 

Without pretending to do justice to Coppersmith's beautiful paper Ref. 0, 
here is a brief proof of his FFT algorithm, Eq.(2d.3), starting from the FFT algorithm 
that is usually found in textbooks. For simplicity, we will only consider the case 
Nb = 4. What follows can be easily generalized to arbitrary A''^ > 1. 

As usual, let 

27r 27r 

a; = exp(?— ) = exp(«— ) . (B.l) 

Define matrices Q,Q^, . . . by 

Q = diag{l,Lj), Q'^ = diag{l,uj'^), ... (B.2) 

Prior to Refs.0, the FFT algorithm might have been stated like this: (see Refs.|0 
and PI) 

F^ = ^{H(E) Is)[Is © (fi^ ® fi' ® fi')](/2 ® F3)P^ , (B.3) 

V2 

Fs = ^{H®h)[h(B{n^®n^)]{l2®F2)P3, (B.4) 

F2 = -^{H ® /2)[/2 © n')]{l2 ® F,)P2 , (B.5) 

F.^^^H. (B.6) 

The matrices P4, P3, P2 are permutation matrices to be specified later. Although we 
could have written just a single equation that combined all four Eqs.(B.3) to (B.6), 
we have chosen not to do this, so as to make explicit the recursive nature of the beast. 
Note 



Q^®Q^®Q = diag{l, u, uj\ uj\ u\ u\ u\ u') = uj^M2)+2Mi)+2'^nio) (g^^) 

The last equation becomes clear when one realizes that 2^n(2) + 2^n(l) + 2°?2(0) 
operating on a state [03, 02, ai, ao) gives the binary expansion of (i(0, 02, ai, a^). Once 
Eq.(B.7) sinks into the old bean, it is only a short step to the realization that: 
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/8©(fi^®fi^®fi) = ci;'^(3)[2 "(2)+2'"(l)+2""(0)] ^ gi[n(3)n(2)02+n(3)n(l)</.3+n(3)n(O)<?i4] ^ (^gj 



h © (^"^ ® fi^) = (CJ^)"(^)[^'"(^)+^""(°)1 = e*["(2)n(l)?i2+n(2)n(0)</.3] ^ j^g^g^ 



/2 © (fi^) = (Lj'')"(l)[20n(0)] ^ gi[n(l)n(O)02] 



(B.IO) 



The final step is to replace all tensor products that contain H by bit operators: 

H^Is = H{3), (B.U) 



h®H®h = H{2) , 



(B.12) 



h®H®l2 = H{1) 



(B.13) 



Is0H = H{0) 



(B.14) 



^^ 



(a^, a^ , a^ , a^) 



K'^'^2'^) 



I(8)P 

w 2^ 3 BR 

(a„ , a_ ,'rra^ ) 
(a , a , a , a ) -^ 

V ' 1 ' 2 ' 3 "^ 
Fig. 11 Permutation matrices P4, P3, P2, Pb_r that arise in the FFT algorithm for 
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If all the permutation matrices of Eqs.(B.3) to (B.6) are combined into a single 
permutation Pbr, then 

Per = {h ® P2){h ® Pz)Pa • (B.15) 

Let a G Bool^ label the columns of a 16 x 16 matrix. As was known previous to 
Ref. 0, the matrices P4,P3,P2 act as shown in Fig.ll. Therefore, their product Pbr 
takes a = (03,02,01,00) to (00,01,02,03); i.e., Prr reverses the bits of o. Pbr is 
not an elementary operation. However, since it is a bit permutation, it can easily be 
expressed as a product of Exchangers, which are elementary. 

Appendix C. Example of Decomposition 
of Controlled Gate into SEO 



Suppose Nb > 3 and consider the following controlled gate: 

g = 5(2)"(^)"(°) . (C.l) 

Replacing n by |[1 — cr^], we get 

g = p(2)3[i--^(i)][i-^Wl . (C.2) 

Expanding the exponent of B{2) in Eq.(C.2), and using identity Eq.(4a.9), we get 

g = AooAoiAi^Aio , (C.3) 

where 

Aoo = i?(2)^ (C.4) 

Aoi = S(2)^"^(°) = P(2)3"(°)S(2)^ , (C.5) 

Au = P(2)3"^W"^(°) = a,(l)"(°)0P(2)^--« = [a,(l)"(°)0P(2)^"(i)]P(2)i , (C.6) 

Aio = P(2)^"^(i) = 5(2)3"«S(2)^ . (C.7) 

The two -8(2)4 and two 5(2)^ factors cancel each other out. We have ordered 
the factors A^ in Eq.(C.3) so that their subscripts b G BooP are in a lazy ordering. 
(Ref. P does this too.) For the above case, lazy ordering is not important. But when 
g has more than two controls, lazy ordering permits some c-nots between adjacent 
Ag's to cancel. The same techniques can be used to decompose a controlled gate with 
more than two controls. 
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Appendix D. Transforming between 

B{n^n) and B{(Jz^l2) bases 

In this appendix, we will show how to transform between the B{n,n) and 
^io'z, h) bases of the set of Ns x Ns diagonal complex matrices. 
Suppose that Nb = 1 and let 



X 



X' 



n(0) " 
n(0) 

1 



Then 



where 



X = MiX' , 



Ml = - 



Note that 



Next suppose that A''^ = 2. Let 



Hi. 



'-H, 
2 



X 



X' 



n(l)n(0) 
n(l)n(0) 
n(l)n(0) 
n(l)n(0) 

1 

a,(l)a,(0) 



n(l) 
n(l) 



1 



n(0) 
n(0) 



1 
a.(0) 



Then 



where 



X = M2X' , 



(D.l) 
(D.2) 

(D.3) 
(D.4) 

(D.5) 



(D.6) 



(D.7) 



{m 
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Mo = ^ 







00 


01 


10 


11 


1 


00 


1 


1 


1 


1 


01 


1 


-1 


1 


-1 


22 
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[D.9) 



Note that 



M^'^ = H2 . 



(D.IO) 



The above results for Nb = 1,2 can be easily generalized to any Nb > 1- One 
can prove using induction that for any A^^, the matrix M7V3 has entries given by 






(D.U) 



where a,b E Zo,Arg_i, a = d{a) and h = d{b). The matrix M^r^ is just the A''^ bit 
Sylvester-Hadamard matrix divided by l/Ns- 

If V and v' are Ns dimensional complex vectors such that 



/T v/ 



v'X = v"X' , 



then 



(D.12) 



implies 



X = Mn,X' 



/ = Ml^v . 



(D.13) 



(D.14) 



Appendix E. Transforming between 

B{n, n) and B{n^ I2) bases 

In this appendix, we will show how to transform between the B{n,n) and 
B{(7z, h) bases of the set of Ns x Ng diagonal complex matrices. 
Suppose that A''^ = 1 and let 



X 



X' 



n(0) 
n(0) 

1 
n(0) 



(E.l) 

(E.2) 



Then 
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X = MiX' , 



(E.3) 



where 



Ml 



Note that 



M^^ -- 
Next suppose that A''^ = 2. Let 



1 -1 
1 



1 1 
1 



X 



X' 



n(l)rl(0) 
rl(l)n(0) 
n(l)n(0) 
n(l)n(0) 

1 

n(0) 

n(l) 

n(l)n(0) 



n(l) 
n(l) 



1 

n(l) 



Then 



where 



X = M2X' , 



M, 
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10 
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-1 


11 











1 



n(0) 
n(0) 



1 
n(0) 



Ml OMi 



(E.4) 



(E.5) 



(E.6) 



(E.7) 



(E.8) 



(E.9) 



Note that 



M,' 



Mf ' (g) Ml 



-1 



(E.IO) 



For any non-negative integers a, b, let a&6 and a Ab denote the same thing as 
they do in the C and C++ computer programming languages. That is, a&6 is the 
result of bitwise and-ing of a and b, and a A 6 is the result of bitwise xor-ing of a and 
b (which is the same as bitwise mod(2) addition of a and b). For any a G Bool'^^ , if 
a = d{d), we define ||a||i = ||a||i = J2^i=o o^- ^^ other words, ||a||i is the number of 
ON bits in a. 
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The above results for A^^ = 1,2 can be easily generalized to any A^^ > 1. One 
can prove using induction that for any A^^, the matrix Mjv^ has entries given by 

iM^s)a,t = (-l)ll"^'ll^(5(a&6, a) , (E.ll) 

where a,b E Zq^Ns-i- 

Appendix F. Central Matrices 

Whose Entries Are Elements of { — 1, 0, 1} 

In this paper, we use bare central matrices whose D matrices have angles 
contained in the interval [0, 90"]. An interesting special case is when these angles are 
all either or 90° exclusively. The entries of such central matrices are elements of 
{-1,0,1}. 

First suppose that the central matrix Q is a single D matrix. Then we can 
express Q as 

Q = expit J2 (t>b(^y®Pb] ' (F-1) 

where the 0^ are all either or 90°. Thus 

Q = [e*f-.(JVB-i)]Ese5^£ = [][e^f-.(^B-i)]P£ ^ (F.2a) 

where 

S = {b\b e Bool^^-\ 0g = ^} . (F.2b) 

Hence, Q can always be expressed as a product of controlled qubit rotations. If 
X = J2b(^s ^b' then Q always acts on more than 2 bits at a time except when (1)X = 1 
or when (2)X = Pb{<y) for some b G Bool and some a G Zq^n^^i such that a ^ Nb — 1- 
In the first case, Q is just a simple qubit rotation. In the second case 

Q = [g*i<x,(JVs-l)]n(a) _ (p,3) 

Next suppose that Q is a direct sum of several D matrices. As explained in 
Section 4(b), such a Q can be obtained by applying a bit permutation matrix to a Q 
that, like the one in Eq.(F.l), is a single D matrix. 

Thus, when Q is a direct sum of one or more D matrices, Q always acts on 
more than 2 bits except when it is a simple qubit rotation or when it has the form 

Q = [e^f-.(7)]n{") ^ (F.4) 
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where a and 7 are distinct elements of Zo^at^-i and h G Bool. Below, we will show 
that Eq.(F.4) implies 

Hence, such a Q can be expressed as a product of a c-not and a controlled phase shift 
with 2 controls. Whenever the optimization of Section 5(c) (Extracting Phases from 
Complex D Matrices) is ON, we will decompose any Q of the form Eq.(F.4) into the 
SEOof Eq.(F.5). 



Note that 



e'2'^!' 



1 
-1 



(F.6) 

[-ir = dtag{-l,l), (F.7) 

'-ir = dzag{l,-l). (F.8) 



Thus, 



e'^''^ = a,{-ir = (-l)V, . (F.9) 

{n changes to n because ax^z = —<^z<^x-) Eq.(F.5) follows directly from Eq.(F.9). 

Note that we could also decompose the Q of Eq.(F.4) using the techniques of 
Section 4 and Appendix C: 

Neither decomposition Eq.(F.5) nor decomposition Eq.(F.lO) uses controlled qubit 
rotations. The SEO of Eq.(F.5) is shorter than the SEO of Eq.(F.lO). Eq.(F.5) uses 
a controlled phase shift whereas Eq.(F.lO) doesn't. Since controlled phase shifts are 
common when the optimization of Section 5(c) is ON, it is natural to use Eq.(F.5) 
when said optimization is ON. 
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FIGURE CAPTIONS: 

FiG.l A CSD binary tree. 

Fig. 2 Pictorial representation of quantum FFT algorithm. For simplicity, we only 
show the Hadamard matrices. White (unshaded) regions inside a matrix repre- 
sent zero entries. 

Fig. 3 Circuit symbols for the 4 different types of c-nots. 

Fig. 4 Four equivalent circuit diagrams for Exchanger. 

Fig. 5 Four equivalent circuit diagrams for Twin-to- twin-er. 

Fig. 6 Circuit symbol for Exchanger. 

Fig. 7 Circuit diagram for Eq.(3b.7). 

Fig. 8 Subset of a CSD tree. 

Fig. 9 Output of Qubiter with input matrix equal to —^Hn^ for Nb = 2, 3, 4. H^g 
is the Nb bit Hadamard matrix defined in Section 2(c). 

Fig. 10 Output of Qubiter with input matrix equal to PbrF^b for A''^ = 2, 3, 4. Fjq^ 
is the Nb bit DFT matrix defined in Section 2(d). 

Fig. 11 Permutation matrices P4^, P;^, P2, Pbr that arise in the FFT algorithm for 

Nb = 4. 
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